Análisis de Series Temporales

IOT Analytics

library(dplyr)
library(tidyverse)
library(kableExtra)
library(DT) # DataTables



library(imputeTS) # Tratamiento de datos faltantes

library(plotly)   # Gráficas
library(ggplot2)  # Gráficas
library(ggfortify)
library(gridExtra)

library(forecast) # Tratamiento de series temporales

load(file="Datos/DF_Energia_GMinutos.RData")
load(file="Datos/DF_Energia_GHoras.RData")
load(file="Datos/DF_Energia_GDiaria.RData")
load(file="Datos/DF_Energia_GMensual.RData")

Objetivo

Proyección del consumo energético a través de series temporales.

Procedimiento

Estudiaremos los datos de consumo energético correspondientes a los años 2007 a 2009 a través de series temporales.

Haremos proyección de futuro para el año 2010 y compararemos los resultados con los datos reales de este último año (serie real).

Series que vamos a estudiar:

Evolución mensual del consumo para los años 2007-2009. Predicción para el año 2010
Evolución diaria del consumo. Años 2007-2009. Predicción para el año 2010
Evolución semanal del consumo. Años 2007-2009 . Predicción para el año 2010

Evolución mensual del consumo. Años 2007-2009

Serie temporal y visualización

Vamos a representar la serie por meses para cada submedidor y tambien para la energía global co granularidad mensual. En este caso, la frecuencia será 12, tenemos una observación por mes para los años 2007,8 y 2009.

seriesMensuales <- Granularidad_meses %>% filter(year != 2006 & year != 2010)
## Creamos un objeto TS
tsSM1_Mensual<- ts(seriesMensuales$Sub_metering_1, frequency=12, start=c(2007,1))
tsSM2_Mensual<- ts(seriesMensuales$Sub_metering_2, frequency=12, start=c(2007,1))
tsSM3_Mensual<- ts(seriesMensuales$Sub_metering_3, frequency=12, start=c(2007,1))
tsSM_GAP_Mensual<- ts(seriesMensuales$Global_active_power, frequency=12, start=c(2007,1))

## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Mensual, ts.colour = 'red', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía de la cocina \n Años 2007-2009")# + scale_x_date(date_labels = '%b%')

autoplot(tsSM2_Mensual, xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía de la lavandería \n Años 2007-2009")

autoplot(tsSM3_Mensual, ts.colour = 'blue', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

autoplot(tsSM_GAP_Mensual, ts.colour = 'green', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía \n Años 2007-2009")

Parece que los datos tienen estacionariedad (movimientos que se repiten cada año), esto ocurre para todas las variables.

Predicción del año 2010 antes de descomponer la serie

Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.

SM1. Cocina

Modelo

fitSM1_Mensual <- tslm(tsSM1_Mensual ~ trend + season) 
summary(fitSM1_Mensual)


Call:
tslm(formula = tsSM1_Mensual ~ trend + season)

Residuals:
     Min       1Q   Median       3Q      Max 
-20522.7  -5190.9    215.7   6726.8  16702.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  66531.9     7294.5   9.121 4.21e-09 ***
trend         -173.1      200.7  -0.862 0.397393    
season2     -19137.6     9635.6  -1.986 0.059061 .  
season3      -1263.5     9641.8  -0.131 0.896879    
season4     -15856.4     9652.3  -1.643 0.114034    
season5      -6134.4     9666.9  -0.635 0.531966    
season6      -9988.0     9685.6  -1.031 0.313160    
season7     -27046.6     9708.4  -2.786 0.010505 *  
season8     -38674.8     9735.4  -3.973 0.000602 ***
season9     -10641.7     9766.3  -1.090 0.287159    
season10    -15040.7     9801.3  -1.535 0.138536    
season11     -6105.3     9840.3  -0.620 0.541071    
season12     -3378.2     9883.2  -0.342 0.735595    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11800 on 23 degrees of freedom
Multiple R-squared:  0.5796,    Adjusted R-squared:  0.3603 
F-statistic: 2.643 on 12 and 23 DF,  p-value: 0.02184

Resultados:

Valor del coeficiente de determinación: 0.5796. Es un valor bajo, el modelo únicamente explica el 58% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 20 días

forecastfitSM1_Mensual <- forecast(fitSM1_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

SM2. Lavadero

Modelo

fitSM2_Mensual <- tslm(tsSM2_Mensual ~ trend + season) 
summary(fitSM2_Mensual)


Call:
tslm(formula = tsSM2_Mensual ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-26325  -4033   1223   6470  16873 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  83389.3     7617.6  10.947 1.35e-10 ***
trend         -898.9      209.6  -4.289 0.000274 ***
season2     -13475.5    10062.5  -1.339 0.193598    
season3       7508.4    10069.0   0.746 0.463406    
season4     -15023.1    10079.9  -1.490 0.149708    
season5      -7727.5    10095.2  -0.765 0.451777    
season6     -11683.0    10114.7  -1.155 0.259931    
season7     -19075.8    10138.6  -1.882 0.072621 .  
season8     -28408.6    10166.7  -2.794 0.010305 *  
season9      -9204.4    10199.1  -0.902 0.376160    
season10      4809.1    10235.6   0.470 0.642890    
season11     -1696.0    10276.3  -0.165 0.870356    
season12     -3191.5    10321.1  -0.309 0.759941    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12320 on 23 degrees of freedom
Multiple R-squared:  0.6541,    Adjusted R-squared:  0.4736 
F-statistic: 3.624 on 12 and 23 DF,  p-value: 0.003891

Resultados:

Valor del coeficiente de determinación: 0.6541. Es un valor bajo, el modelo únicamente explica el 65% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM2_Mensual <- forecast(fitSM2_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

SM3. AC y termo

Modelo

fitSM3_Mensual <- tslm(tsSM3_Mensual ~ trend + season) 
summary(fitSM3_Mensual)


Call:
tslm(formula = tsSM3_Mensual ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-87112 -18443   4781  17647  82107 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   299382      23189  12.911 5.06e-12 ***
trend           1866        638   2.926  0.00761 ** 
season2       -50805      30631  -1.659  0.11077    
season3       -29034      30651  -0.947  0.35336    
season4       -64190      30684  -2.092  0.04768 *  
season5       -54165      30730  -1.763  0.09125 .  
season6       -85197      30790  -2.767  0.01097 *  
season7      -144850      30863  -4.693 1.00e-04 ***
season8      -171144      30948  -5.530 1.27e-05 ***
season9       -69624      31047  -2.243  0.03486 *  
season10      -52696      31158  -1.691  0.10430    
season11      -37281      31282  -1.192  0.24550    
season12        7092      31418   0.226  0.82341    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 37510 on 23 degrees of freedom
Multiple R-squared:  0.7577,    Adjusted R-squared:  0.6312 
F-statistic: 5.992 on 12 and 23 DF,  p-value: 0.0001238

Resultados:

Valor del coeficiente de determinación: 0.7577. Es un valor que podría ser aceptable, el modelo únicamente explica el 76% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM3_Mensual <- forecast(fitSM3_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

Comparación de los coeficientes y resultados de cada modelo

DatosRealesSeriesMensuales <- Granularidad_meses %>% filter(year == 2010 )

RMSE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
R2_SM1_GranMensual<-0.5796
MASE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_1)[12]

RMSE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
R2_SM2_GranMensual<-0.6541
MASE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_2)[12]



RMSE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
R2_SM3_GranMensual<-0.7577
MASE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_3)[12]

ResEvolMensual<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensual,RMSE_SM2_GranMensual,RMSE_SM3_GranMensual),
     # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
      R2 = c(R2_SM1_GranMensual,R2_SM2_GranMensual,R2_SM3_GranMensual)
)

ResEvolMensual %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE	R2
Cocina	11186.74	0.5796
Lavadero	10939.60	0.6541
AC y TermoE	57300.05	0.7577

# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.

El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.

Forecasting descomponiendo la serie (así ok)

Descomponer la serie temporal: identificar y calcular las componentes tendencia, estacionalidad y ruido.

Tendencia de una serie temporal: Consiste en la evolución a largo plazo de la serie. En nuestro caso, puede haber una tendencia al aumento del consumo energético o una tendencia decreciente.
Estacionalidad: Movimientos oscilatorios alrededor de la tendencia que se repiten de forma periódica y con amplitud inferior al año
Componente irregular: es aleatoria y no se puede predecir (ruido)

Descomposición de la serie y visualización

Cocina: Submetering 1

Componentes_SM1_MensualSTL<- tsSM1_Mensual %>%stl( t.window=12,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_MensualSTL, ts.colour = 'red', main="Descomposición de componentes")

Se aprecia una gran influencia de la componente estacional: oscilaciones que se repiten cada año, por ejemplo, en el mes de agosto, el consumo energético de la cocina decrece consuderablemente y es mucho menor que el resto de meses. Por el contrario, es en el mes de Junio cuando se registra el mayor consumo energético de los aparatos submeterizados de la cocina.

Tambien se aprecia una tendencia creciente de consumo, cada año se consume más.

Remainder no es un muelle, gráficamente, parece que la componente irregular no es del todo aleatoria.

plot.ts(tsSM1_Mensual, col = 'gray',
         xlab = "Mes del año", 
        ylab = "Energía (Varios-hora)", 
        main = "Evolución mensual del consumo energético de la cocina \n y la tendencia del consumo")
lines(Componentes_SM1_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Como ya habíamos comentado, existe una tendencia creciente del consumo energético de la cocina.

Lavadero: Submetering 2

Componentes_SM2_MensualSTL<- tsSM2_Mensual %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_MensualSTL, ts.colour = 'blue')

Por el contrario que en la cocina, se observa como existe una considerable tendencia decreciente del consumo en el lavadero.

Tambien se observa claramente la estacionalidad, donde podemos apreciar que en los meses de verano el consumo energético es mucho menor que el resto de meses del año.

Parece que la componente irregular podría ser aleatoria, sin embargo, esto es únicamente a modo gráfico y por tanto deberíamos comprobarlo con un test de hipótesis.

plot.ts(tsSM2_Mensual, col = 'gray', main="Evolución mensual del consumo energético del lavadero")
lines(Componentes_SM2_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2
      )

Termo eléctrico y aire acondicionado: Submetering 3

tsSM3_Mensual %>%
  stl( t.window=12,
      s.window="periodic", 
      robust=TRUE)  %>% summary()

 Call:
 stl(x = ., s.window = "periodic", t.window = 12, robust = TRUE)

 Time.series components:
    seasonal              trend            remainder        
 Min.   :-107286.17   Min.   :236737.3   Min.   :-98598.49  
 1st Qu.: -12301.07   1st Qu.:257514.6   1st Qu.:-14547.16  
 Median :   8318.49   Median :277783.5   Median :  -493.40  
 Mean   :      0.00   Mean   :274647.8   Mean   : -3392.79  
 3rd Qu.:  27345.56   3rd Qu.:287394.0   3rd Qu.:  9934.48  
 Max.   :  72770.26   Max.   :312023.9   Max.   : 91583.87  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     39647        29879     24482         82030
   %  48.3         36.4      29.8         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8138  0.9452  0.8280  0.9944  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 361 12 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 37 2 2
 $ inner: int 1
 $ outer: int 15

Componentes_SM3_MensualSTL<- tsSM3_Mensual %>%
  stl( t.window=12, s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_MensualSTL, ts.colour = 'blue')

Tendencia creciente al consumo energético. Se observa fuerte componente estacional, donde el consumo energético baja considerablemente en los meses de verano, y alcanza niveles máximos en los últimos meses del año (Noviembre y Diciembre)

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

plot.ts(tsSM3_Mensual, col = 'gray',
        main="Evolución del consumo energético mensual del aire acondicionado y termo eléctrico")
lines(Componentes_SM3_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Vamos a aplicar suavizado exponencial para realizar un pronóstico del comportamiento de consumo energético. Veremos esta predicción para datos:

Sin estacionalidad
Con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes

Evolución mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina

tsSM1_Mensual_detrended <- tsSM1_Mensual - Componentes_SM1_MensualSTL$time.series[,1]
# Elimino la componente estacional
plot.ts(tsSM1_Mensual_detrended, 
        ylab = "Energía (Vatios-hora)", xlab="Mes del año", #ylim=c(0,3000),
        main = 'Evolución del consumo energético de la cocina sin estacionalidad',
        cex.main = 0.85, col="gray")

Vamos a estudiar como se comporta la componente irregular.

par(mfrow = c(1,2))
plot(Componentes_SM1_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_MensualSTL$time.series[,3])
qqline(Componentes_SM1_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM1_MensualSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM1_MensualSTL$time.series[, 3]
W = 0.75719, p-value = 2.554e-06

Los resíduos no son solo ruido. Aún hay estacionalidad.

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Mensual_detrended ,s.window = "periodic",robust = TRUE),colour = 'red', main="Descomposición de la serie desestacionalizada")

Podríamos a sumir que se ha eliminado la estacionalidad, ya que los resíduos se comportan de manera irregular.

Suavizado exponencial de HoltWinters

Primera iteracción: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicción

tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
#   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)

      xhat           level           trend            season     
 Min.   : 7683   Min.   :13713   Min.   :-376.9   Min.   :-7702  
 1st Qu.:45514   1st Qu.:50966   1st Qu.:-376.9   1st Qu.:-5695  
 Median :49998   Median :52578   Median :-376.9   Median :-1479  
 Mean   :49432   Mean   :49809   Mean   :-376.9   Mean   :    0  
 3rd Qu.:57209   3rd Qu.:58150   3rd Qu.:-376.9   3rd Qu.: 1006  
 Max.   :88889   Max.   :66538   Max.   :-376.9   Max.   :25965

# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronóstico de consumo será en el rango: 7683-88889 vatios-hora de energía consumida.

Segunda iteracción: Añadiendo \(\beta\) podré averiguar la tendencia del consumo para las predicciones

tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)

      xhat            level           trend              season     
 Min.   :-19652   Min.   :13713   Min.   :-36810.3   Min.   :-7702  
 1st Qu.: 44754   1st Qu.:50966   1st Qu.: -3854.5   1st Qu.:-5695  
 Median : 50051   Median :52578   Median :  -188.6   Median :-1479  
 Mean   : 50152   Mean   :49809   Mean   :   343.0   Mean   :    0  
 3rd Qu.: 69761   3rd Qu.:58150   3rd Qu.:  5871.7   3rd Qu.: 1006  
 Max.   :101263   Max.   :66538   Max.   : 42426.8   Max.   :25965

# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia al consumo decreciente, es decir, se pronosticará un descenso del consumo par meses siguientes.

Tercera iteracción: Si añado \(\gamma\) podré ver la estacionalidad.

tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
 # gamma=TRUE,
   alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)

      xhat            level           trend              season     
 Min.   :-19652   Min.   :13713   Min.   :-36810.3   Min.   :-7702  
 1st Qu.: 44754   1st Qu.:50966   1st Qu.: -3854.5   1st Qu.:-5695  
 Median : 50051   Median :52578   Median :  -188.6   Median :-1479  
 Mean   : 50152   Mean   :49809   Mean   :   343.0   Mean   :    0  
 3rd Qu.: 69761   3rd Qu.:58150   3rd Qu.:  5871.7   3rd Qu.: 1006  
 Max.   :101263   Max.   :66538   Max.   : 42426.8   Max.   :25965

# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

# Rstudio optimiza los coeficientes
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)

      xhat           level           trend            season       
 Min.   :41445   Min.   :46024   Min.   :-376.9   Min.   :-7930.2  
 1st Qu.:44274   1st Qu.:48191   1st Qu.:-376.9   1st Qu.:-5667.3  
 Median :47013   Median :50358   Median :-376.9   Median :-1903.9  
 Mean   :49704   Mean   :50358   Mean   :-376.9   Mean   : -276.8  
 3rd Qu.:51922   3rd Qu.:52525   3rd Qu.:-376.9   3rd Qu.: 1093.8  
 Max.   :76135   Max.   :54692   Max.   :-376.9   Max.   :25965.5

# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_for <- forecast(tsSM1_HW_Mensual, h=12)
plot(tsSM1_HW_Mensual_for ,
     ylab= "Vatios-Hora",
     xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forC <- forecast(tsSM1_HW_Mensual, h=12
                                  , 
                                  level=c(10,25)
                                  )
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha",
     main="Pronóstico del consumo de la cocina. Año 2010 \n (sin estacionalidad)", start(2010))

Lavadero

tsSM2_Mensual_detrended <- tsSM2_Mensual - Componentes_SM2_MensualSTL$time.series[,2]
plot.ts(tsSM2_Mensual_detrended, 
        ylab = "Energía (Vatios-hora)", col="gray",
        main = 'Variabilidad del consumo energético del lavadero \n para datos sin estacionalidad',
        cex.main = 0.85)

Vamos a estudiar como se comporta la componente irregular.

par(mfrow = c(1,2))
plot(Componentes_SM2_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_MensualSTL$time.series[,3])
qqline(Componentes_SM2_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM2_MensualSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM2_MensualSTL$time.series[, 3]
W = 0.85783, p-value = 0.0002898

Los resíduos no siguen una distribución normal. La descomposición no es buena, por tanto, tenemos que vovler a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Mensual_detrended ,s.window = "periodic",
         robust = TRUE,
         t.window = 12))

Ahora si que parece que los resíduos se comportan de manera aleatoria. Comportamiento de consumo: Hay una tendencia decreciente al consumo los 5 primeros meses del año 2007, sin embargo, posteriormente este comportamiento cambia y existe una tendencia al consumo hasta principios del año 2008, donde vuelve a cambiar la tendencia y pasa a haber un comportamiento decreciente hasya finales del año 2008. Al parecer, desde final del año 2008 se ha observado una tendencia creciente del consumo energético del lavadero.

Suavizado exponencial de HoltWinters

Primera iteracción: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicción

tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
#   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)

      xhat              level            trend           season      
 Min.   :-37423.2   Min.   :-21787   Min.   :237.4   Min.   :-22162  
 1st Qu.: -9524.1   1st Qu.: -7651   1st Qu.:237.4   1st Qu.: -8854  
 Median : -1669.0   Median :  1753   Median :237.4   Median :  3244  
 Mean   :  -131.6   Mean   :  -369   Mean   :237.4   Mean   :     0  
 3rd Qu.: 11559.8   3rd Qu.:  3709   3rd Qu.:237.4   3rd Qu.:  8595  
 Max.   : 27452.5   Max.   : 32063   Max.   :237.4   Max.   : 17041

# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronóstico de consumo será en el rango: 0-27452.5 vatios-hora de energía consumida, ya que no puede haber un consumo energético negativo.

Segunda iteracción: Añadiendo \(\beta\) podré averiguar la tendencia del consumo para las predicciones

tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)

      xhat              level            trend              season      
 Min.   :-62015.7   Min.   :-21787   Min.   :-27388.9   Min.   :-22162  
 1st Qu.:-13466.1   1st Qu.: -7651   1st Qu.: -6777.2   1st Qu.: -8854  
 Median :  2156.0   Median :  1753   Median :  -158.7   Median :  3244  
 Mean   :  -144.9   Mean   :  -369   Mean   :   224.2   Mean   :     0  
 3rd Qu.: 16152.5   3rd Qu.:  3709   3rd Qu.:  7242.4   3rd Qu.:  8595  
 Max.   : 73833.0   Max.   : 32063   Max.   : 53850.3   Max.   : 17041

# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia al consumo cambiante como ya hemos comentado antes

Tercera iteracción: Si añado \(\gamma\) podré ver la estacionalidad.

tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)

      xhat              level            trend              season      
 Min.   :-62015.7   Min.   :-21787   Min.   :-27388.9   Min.   :-22162  
 1st Qu.:-13466.1   1st Qu.: -7651   1st Qu.: -6777.2   1st Qu.: -8854  
 Median :  2156.0   Median :  1753   Median :  -158.7   Median :  3244  
 Mean   :  -144.9   Mean   :  -369   Mean   :   224.2   Mean   :     0  
 3rd Qu.: 16152.5   3rd Qu.:  3709   3rd Qu.:  7242.4   3rd Qu.:  8595  
 Max.   : 73833.0   Max.   : 32063   Max.   : 53850.3   Max.   : 17041

# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

# Rstudio optimiza los coeficientes
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
 #  alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)

      xhat               level             trend           season        
 Min.   :-28433.55   Min.   :-1983.2   Min.   :234.8   Min.   :-31197.1  
 1st Qu.: -7760.62   1st Qu.: -612.7   1st Qu.:236.3   1st Qu.: -9057.5  
 Median :  2416.46   Median :  755.1   Median :237.5   Median :  1068.6  
 Mean   :    86.66   Mean   :  751.1   Mean   :237.2   Mean   :  -901.7  
 3rd Qu.:  9550.00   3rd Qu.: 2117.5   3rd Qu.:238.2   3rd Qu.:  8271.5  
 Max.   : 18504.81   Max.   : 3472.5   Max.   :239.1   Max.   : 17041.2

# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_for <- forecast(tsSM2_HW_Mensual, h=12)
plot(tsSM2_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forC <- forecast(tsSM2_HW_Mensual, h=12, level=95)
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético del lavadero. Año 2010", start(2010))

AC y termo eléctrico

tsSM3_Mensual_detrended <- tsSM3_Mensual - Componentes_SM3_MensualSTL$time.series[,2]
plot.ts(tsSM3_Mensual_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Evolución el consumo energético del aire acondicionado y termo ele´ctrico \n Datos sin estacionalidad',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_MensualSTL$time.series[,3])
qqline(Componentes_SM3_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM3_MensualSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM3_MensualSTL$time.series[, 3]
W = 0.87556, p-value = 0.0007821

Los resíduos no siguen una distribución normal. La descomposición no es buena.

Hay que volver a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Mensual_detrended ,s.window = "periodic",robust = TRUE))

Parece que la componente irregular es aleatoria. La tendencia del consumo energético es creciente, aunque durante el año 2008 hubo períodos con tendencia decreciente, al igual que los 5 primeros meses del año 2007.

Suavizado exponencial de HoltWinters

Primera iteracción: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicción

tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , seasonal = "additive",
#   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM3_HW_Mensual$fitted)

      xhat             level             trend          season      
 Min.   :-209442   Min.   :-181400   Min.   :-635   Min.   :-92079  
 1st Qu.:  -9333   1st Qu.:  -8991   1st Qu.:-635   1st Qu.:-23534  
 Median :   9306   Median :   4859   Median :-635   Median :   699  
 Mean   :  -5753   Mean   :  -5118   Mean   :-635   Mean   :     0  
 3rd Qu.:  31822   3rd Qu.:  20231   3rd Qu.:-635   3rd Qu.: 17164  
 Max.   :  80208   Max.   :  45472   Max.   :-635   Max.   : 86727

# optimiza
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronóstico de consumo será en el rango: 0-80208 vatios-hora de energía consumida, ya que no puede haber un consumo energético negativo.

Segunda iteracción: Añadiendo \(\beta\) podré averiguar la tendencia del consumo para las predicciones

tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM3_HW_Mensual$fitted)

      xhat             level             trend               season      
 Min.   :-422399   Min.   :-181400   Min.   :-213592.3   Min.   :-92079  
 1st Qu.: -27389   1st Qu.:  -8991   1st Qu.: -26879.7   1st Qu.:-23534  
 Median :  11408   Median :   4859   Median :  -1375.6   Median :   699  
 Mean   :  -5558   Mean   :  -5118   Mean   :   -439.2   Mean   :     0  
 3rd Qu.:  53303   3rd Qu.:  20231   3rd Qu.:  15913.9   3rd Qu.: 17164  
 Max.   : 245407   Max.   :  45472   Max.   : 213552.8   Max.   : 86727

# optimiza
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia al consumo cambiante como ya hemos comentado antes.

Tercera iteracción: Si añado \(\gamma\) podré ver la estacionalidad.

tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)

summary(tsSM3_HW_Mensual$fitted)

      xhat             level             trend               season      
 Min.   :-422399   Min.   :-181400   Min.   :-213592.3   Min.   :-92079  
 1st Qu.: -27389   1st Qu.:  -8991   1st Qu.: -26879.7   1st Qu.:-23534  
 Median :  11408   Median :   4859   Median :  -1375.6   Median :   699  
 Mean   :  -5558   Mean   :  -5118   Mean   :   -439.2   Mean   :     0  
 3rd Qu.:  53303   3rd Qu.:  20231   3rd Qu.:  15913.9   3rd Qu.: 17164  
 Max.   : 245407   Max.   :  45472   Max.   : 213552.8   Max.   : 86727

# optimiza
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El ajuste podría ser mejorable, ya que hay gran influencia de los valores extremos.

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_for <- forecast(tsSM3_HW_Mensual, h=12)
plot(tsSM3_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forC <- forecast(tsSM3_HW_Mensual, h=12, level=95)
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción el consumo energético del lavadero y termo eléctrico para el año 2010", start(2010),
     ylim=c(0,80000))

Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad

En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.

RMSE_SM1_GranMensualFor<-accuracy(tsSM1_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualFor<-accuracy(tsSM2_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualFor<-accuracy(tsSM3_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualFor,
               RMSE_SM2_GranMensualFor,
               RMSE_SM3_GranMensualFor)
)
ResEvolMensual2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE
Cocina	14278.84
Lavadero	44294.27
AC y TermoE	325349.90

Evolución mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.

Cocina

tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)

## Volvemos a descomponer la serie
plot(stl(tsSM1_Mensual_AdjustedII,s.window = "periodic",robust = TRUE,t.window = 12))

Tendencia creciente del consumo.

Suavizado exponencial de HoltWinters

Primera interacción: para observar el rango de la predicción.

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, alpha=TRUE,
                                #  beta=TRUE, gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

- Segunda interacción: añado \(\beta\) para observar la tendencia del consumo:

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, alpha=TRUE,
                                  beta=TRUE # , gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

- Tercera iteracción: Añado \(\gamma\)

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII
                                  , 
                                  alpha=TRUE,
                                  beta=TRUE  , gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, 
                                 #  alpha=TRUE,
                                   beta=TRUE  ,
                                  gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Si dejamos que la propia función optimice los valores de beta y gamma, la predicción se ajusta mejor y por tanto, obtendremos mejores resultados.

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=12, level=95)
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha", main="Pronóstico para el año 2010 (con estacionalidad)", start(2010))

Lavadero

tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)

## Volvemos a descomponer la serie
plot(stl(tsSM2_Mensual_AdjustedII,s.window = "periodic",
         robust = TRUE,t.window = 12))

LA tendencia de consumo es claramente creciente, y parece que la componente irregular se ha eliminado.

Suavizado exponencial de HoltWinters

Primera interacción: componente alpha para ver el rango de predicciones

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                  # beta=TRUE, 
                                  # gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM2_HW_MensualII$fitted)

      xhat           level           trend           season      
 Min.   :10670   Min.   :27252   Min.   :-1121   Min.   :-22206  
 1st Qu.:41144   1st Qu.:44247   1st Qu.:-1121   1st Qu.: -8213  
 Median :56659   Median :56542   Median :-1121   Median :  3157  
 Mean   :52598   Mean   :53719   Mean   :-1121   Mean   :     0  
 3rd Qu.:64844   3rd Qu.:62657   3rd Qu.:-1121   3rd Qu.:  8459  
 Max.   :80481   Max.   :80053   Max.   :-1121   Max.   : 17218

plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que la predicción del consumo energético del lavadero estará entre 10670 y 80481 vatios-hora de energía.

Segunda interaccióm: añado beta para ver la tendencia del consumo:

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                   beta=TRUE, 
                                  # gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM2_HW_MensualII$fitted)

      xhat            level           trend              season      
 Min.   :-12375   Min.   :27252   Min.   :-27249.7   Min.   :-22206  
 1st Qu.: 39497   1st Qu.:44247   1st Qu.: -7920.3   1st Qu.: -8213  
 Median : 56536   Median :56542   Median : -2034.5   Median :  3157  
 Mean   : 52768   Mean   :53719   Mean   :  -951.5   Mean   :     0  
 3rd Qu.: 66700   3rd Qu.:62657   3rd Qu.:  6611.9   3rd Qu.:  8459  
 Max.   :121207   Max.   :80053   Max.   : 52800.2   Max.   : 17218

plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Como ya habíamos avanzado, tendencia creciente

Tercera interacción: añadimos gamma

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                   beta=TRUE, 
                                   gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM2_HW_MensualII$fitted)

      xhat            level           trend              season      
 Min.   :-12375   Min.   :27252   Min.   :-27249.7   Min.   :-22206  
 1st Qu.: 39497   1st Qu.:44247   1st Qu.: -7920.3   1st Qu.: -8213  
 Median : 56536   Median :56542   Median : -2034.5   Median :  3157  
 Mean   : 52768   Mean   :53719   Mean   :  -951.5   Mean   :     0  
 3rd Qu.: 66700   3rd Qu.:62657   3rd Qu.:  6611.9   3rd Qu.:  8459  
 Max.   :121207   Max.   :80053   Max.   : 52800.2   Max.   : 17218

plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                   beta=TRUE, 
                                   gamma=TRUE #,
                                 #  alpha = TRUE
                                  )

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=12)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=12, level=95)
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico

tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)

## Volvemos a descomponer la serie
plot(stl(tsSM3_Mensual_AdjustedII,s.window = "periodic",robust = TRUE, 
         t.window = 12))

Tendencia creciente del consumo de este submedidor.

Suavizado exponencial de HoltWinters

Primera interacción: alpha para ver el margen de predicción

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII,
                                #  beta=TRUE, 
                                #  gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM3_HW_MensualII$fitted)

      xhat            level            trend          season       
 Min.   : 80833   Min.   :111338   Min.   :2186   Min.   :-104479  
 1st Qu.:262195   1st Qu.:268973   1st Qu.:2186   1st Qu.: -24576  
 Median :295777   Median :280868   Median :2186   Median :   6764  
 Mean   :281117   Mean   :278931   Mean   :2186   Mean   :      0  
 3rd Qu.:319558   3rd Qu.:307632   3rd Qu.:2186   3rd Qu.:  21778  
 Max.   :392649   Max.   :331536   Max.   :2186   Max.   :  86499

plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Rrango de predicciónes: 80833-392694 vatios hora.

Segunda interacción: añado beta para ver la tendencia del consumo

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII,
                                  beta=TRUE, 
                                #  gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM3_HW_MensualII$fitted)

      xhat             level            trend             season       
 Min.   :-141552   Min.   :111338   Min.   :-220198   Min.   :-104479  
 1st Qu.: 268345   1st Qu.:268973   1st Qu.: -28339   1st Qu.: -24576  
 Median : 293751   Median :280868   Median :   1066   Median :   6764  
 Mean   : 281347   Mean   :278931   Mean   :   2417   Mean   :      0  
 3rd Qu.: 352536   3rd Qu.:307632   3rd Qu.:  23657   3rd Qu.:  21778  
 Max.   : 524064   Max.   :331536   Max.   : 207838   Max.   :  86499

plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que los valores extremos influirán en las predicciones

Tercera interacción

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII,
                                  beta=TRUE, 
                                  gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM3_HW_MensualII$fitted)

      xhat             level            trend             season       
 Min.   :-141552   Min.   :111338   Min.   :-220198   Min.   :-104479  
 1st Qu.: 268345   1st Qu.:268973   1st Qu.: -28339   1st Qu.: -24576  
 Median : 293751   Median :280868   Median :   1066   Median :   6764  
 Mean   : 281347   Mean   :278931   Mean   :   2417   Mean   :      0  
 3rd Qu.: 352536   3rd Qu.:307632   3rd Qu.:  23657   3rd Qu.:  21778  
 Max.   : 524064   Max.   :331536   Max.   : 207838   Max.   :  86499

plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=12)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=12, level=95)
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético del aire acondicionado \n  y termo eléctrico  para el año 2010", start(2010), 
     ylim=c(0,400000))

Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad

RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualForII,
               RMSE_SM2_GranMensualForII,
               RMSE_SM3_GranMensualForII)
   
    
)

ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE
Cocina	16362.16
Lavadero	27223.90
AC y TermoE	83037.90

Comparación de resultados. Predicción antes y despues de eliminar estacionalidad

Comparacion<-cbind.data.frame(
  ResEvolMensual2_SinEstacionalidadConTendencia,
  ResEvolMensual2SinTendenciaConEstacionalidad[,2])

colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE Con Estacionalidad	RMSE Sin estacionalidad
Cocina	14278.84	16362.16
Lavadero	44294.27	27223.90
AC y TermoE	325349.90	83037.90

Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.

Comparación de resultados. Predicción antes y despues de eliminar la tendencia (descomponer)

En las gráficas de predicción de suaviazado exponencial, hemos observado que el ajuste cuando eliminamos la estacionalidad, no mantiene la tendencia de consumo energético de la serie. Sin embargo, al mantener la estacionalidad se puede apreciar como las predicciones ajustan la tendencia de la serie original

# Predicción antes y despues de eliminar estacionalidad:
Comparacion

   Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
1      Cocina                14278.84                16362.16
2    Lavadero                44294.27                27223.90
3 AC y TermoE               325349.90                83037.90

# Sin descomponer


ComparacionII<-cbind.data.frame(ResEvolMensual,
                              Comparacion[-1]
                              )
colnames(ComparacionII) <- c("Submedidor","RMSE","R^2","Con estacionalidad","Sin estacionalidad")
ComparacionII %>% kable(booktabs=TRUE) %>% 
  kable_styling(latex_options = "striped") %>% 
  add_header_above( c("", "Serie sin descomponer" = 2, "Descomponiendo la serie" = 2)) %>% 
  add_header_above( c("", "RMSE" = 4))

	RMSE
	Serie sin descomponer		Descomponiendo la serie
Submedidor	RMSE	R^2	Con estacionalidad	Sin estacionalidad
Cocina	11186.74	0.5796	14278.84	16362.16
Lavadero	10939.60	0.6541	44294.27	27223.90
AC y TermoE	57300.05	0.7577	325349.90	83037.90

Por tanto, nos quedaremos con los modelos correspondientes a la serie sin estacionalidad, que son los que mejores resultados nos ha dado.

Pronóstico de consumo para los próximos 12 meses

Cocina

# Definición de serie temporal
seriesMensualesCompletas <- Granularidad_meses 
## Creamos un objeto TS
tsSM1_Mensual<- ts(seriesMensualesCompletas$Sub_metering_1, frequency=12, start=c(2006,12))
tsSM2_Mensual<- ts(seriesMensualesCompletas$Sub_metering_2, frequency=12, start=c(2006,12))
tsSM3_Mensual<- ts(seriesMensualesCompletas$Sub_metering_3, frequency=12, start=c(2006,12))
tsSM_GAP_Mensual<- ts(seriesMensualesCompletas$Global_active_power, frequency=12, start=c(2006,12))

# Descomposición de la serie
Componentes_SM1_MensualSTL<-  tsSM1_Mensual %>%
  stl( t.window=12,s.window="periodic",  robust=TRUE)
# Eliminación de la estacionalidad
tsSM1_Mensual_detrended <- tsSM1_Mensual - Componentes_SM1_MensualSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , 
                                seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalCocinaAnual <- forecast(tsSM1_HW_Mensual, h=13)
## Plot only the forecasted area
plot(PrediccionFinalCocinaAnual, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético \n  de la cocina para los próximos 12 meses", start(2010,12),
     ylim=c(0,80000))

Las predicciones son las siguientes:

PrediccionFinalCocinaAnual[["mean"]]

          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2010                                                                        
2011 41945.83 39455.73 37523.93 44392.00 45300.17 34904.71 39853.79 39158.80
          Sep      Oct      Nov      Dec
2010                            65059.25
2011 42425.45 39842.72 39116.19 62084.64

Lavadero

# Descomposición de la serie
Componentes_SM2_MensualSTL<-  tsSM2_Mensual %>%
  stl( t.window=12,s.window="periodic",  robust=TRUE)
# Eliminación de la estacionalidad
tsSM2_Mensual_detrended <- tsSM2_Mensual - Componentes_SM2_MensualSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , 
                                seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalLavaderoAnual <- forecast(tsSM2_HW_Mensual, h=13)
## Plot only the forecasted area
plot(PrediccionFinalLavaderoAnual, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético \n  del lavadero para los próximos 12 meses", start(2010,12),
     ylim=c(0,80000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoAnual[["mean"]]

          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2010                                                                        
2011 34865.05 38962.01 33876.42 30611.39 37750.86 30764.14 35867.64 33948.36
          Sep      Oct      Nov      Dec
2010                            56335.97
2011 35359.50 31508.32 31995.30 45813.51

Aire acondicionado y termo eléctrico

# Descomposición de la serie
Componentes_SM3_MensualSTL<-  tsSM3_Mensual %>%
  stl( t.window=12,s.window="periodic",  robust=TRUE)
# Eliminación de la estacionalidad
tsSM3_Mensual_detrended <- tsSM3_Mensual - Componentes_SM3_MensualSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , 
                                seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalAcTAnual <- forecast(tsSM3_HW_Mensual, h=13)
## Plot only the forecasted area
plot(PrediccionFinalAcTAnual, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético \n  del aire y termo eléctricopara los próximos 12 meses", start(2010,12))

Las predicciones son las siguientes:

PrediccionFinalLavaderoAnual[["mean"]]

          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2010                                                                        
2011 34865.05 38962.01 33876.42 30611.39 37750.86 30764.14 35867.64 33948.36
          Sep      Oct      Nov      Dec
2010                            56335.97
2011 35359.50 31508.32 31995.30 45813.51

Evolución semanal del consumo. Años 2007-2009

Serie temporal y visualización

Vamos a representar la serie por semanas. Es decir, para cada submedidor vamos a representar el consumo total semanal. En este caso, la frecuencia será de 7 días por año, tenemos una observación para cada día de la semana para los años 2007,8 y 2009.

Granularidad_semanas <- Granularidad_horas %>% group_by(year, weekday) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )
seriesSemanales <- rbind.data.frame(
   filter(Granularidad_semanas, year== 2007  ),
   filter(Granularidad_semanas, year== 2008  ),
   filter(Granularidad_semanas, year== 2009  )
)

## Creamos un objeto TS
 tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1,frequency=7, start=c(2007,1))
 tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2,frequency=7, start=c(2007,1))
 tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,frequency=7, start=c(2007,1))

## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la cocina \n Años 2007-2009")

autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la lavandería \n Años 2007-2009")

autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

Forecasting descomponiendo la serie

Descomposición de la serie y visualización

Cocina: Submetering 1

Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%
  stl( t.window=7, s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')

Se aprecia componente estacional: oscilaciones que se repiten cada año, es decir, cada 7 días. Vemos como la tendencia de consumo es creciente para los años 2007 y 2009 pero sin embargo, en el año 2008 hubo tendencia de consumo decreciente.

plot.ts(tsSM1_Semanal, col = 'gray',main="Evolución del consumo energético de cada día de la semana de la cocina \n con la tendencia")
lines(Componentes_SM1_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Lavadero: Submetering 2

Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( t.window=7,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')

Se observa una tendencia al consumo que es decreciente respecto de los últimos años, aunque a partir del año 2009 se puede apreciar cierta tendencia creciente del consumo energético del lavadero.

Tambien se obserca componente estacional, donde hay un pico de menor consumo energético los jueves y domingos, y los sábados son lso días donde hay un mayor registro del consumo energético del lavadero.

plot.ts(tsSM2_Semanal, col = 'gray' , main = 'Evolución del consumo energético del lavadero \n por día de la semana. Años 2007-2009')
lines(Componentes_SM2_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Aire acondicionado y termo eléctrico: Submetering 3

Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( t.window=7, 
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')

Se observa una fuerte componente estacional, con patrones que se repiten para cada año los mismos días de la semana. Por ejemplo, los Martes y Viernes son los días donde el consumo energético es mayor y los sábados y Miércoles son los días donde el consumo energético es menor.

La tendencia general del consumo energético es creciente.

plot.ts(tsSM3_Semanal, col = 'gray' , main = 'Evolución del consumo energético del AC y termo \n  según día de la semana. Años 2007-2009',)
lines(Componentes_SM3_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Tendencia creciente del consumo energético del aire acondicionado y termo eléctrico

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Evolución semanal del consumo (sin estacionalidad). Años 2007-2009

Cocina

tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,1]
plot.ts(tsSM1_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Evolución semanal del consumo de la cocina. Años 2007-2009',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.62627, p-value = 3.78e-06

Los resíduos no son solo ruido. Aún hay estacionalidad. Volvemos a descomponer la serie desestacionalizada una vez para tratar de eliminar la componente estacional.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,
         s.window = "periodic",robust = TRUE,t.window = 7))

Parece que la componente irregular no sigue ningún patrón determinado.

Suavizado exponencial de HoltWinters

Primera iteracción: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicción

tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                                alpha=TRUE
                               # beta=TRUE, 
                               # gamma=TRUE
                                )
summary(tsSM1_HW_Semanal$fitted)

      xhat            level            trend            season        
 Min.   : 47685   Min.   : 67975   Min.   :-900.3   Min.   :-19389.3  
 1st Qu.: 81771   1st Qu.: 80568   1st Qu.:-900.3   1st Qu.: -1906.8  
 Median : 85529   Median : 87483   Median :-900.3   Median :   336.8  
 Mean   : 83958   Mean   : 84859   Mean   :-900.3   Mean   :     0.0  
 3rd Qu.: 92100   3rd Qu.: 89084   3rd Qu.:-900.3   3rd Qu.:  6619.6  
 Max.   :105496   Max.   :101417   Max.   :-900.3   Max.   :  9648.0

# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronóstico de consumo de la cocina variará con un pronóstico de 47685 vatios hora como mínimo y 105496 como máximo.

Segunda iteracción: Añadiendo \(\beta\) podré averiguar la tendencia del consumo para las predicciones

tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                                alpha=TRUE,
                                beta=TRUE
                               # gamma=TRUE
                                )
summary(tsSM1_HW_Semanal$fitted)

      xhat            level            trend              season        
 Min.   : 36880   Min.   : 67975   Min.   :-32119.9   Min.   :-19389.3  
 1st Qu.: 73545   1st Qu.: 80568   1st Qu.: -6256.7   1st Qu.: -1906.8  
 Median : 86848   Median : 87483   Median :  -534.6   Median :   336.8  
 Mean   : 83905   Mean   : 84859   Mean   :  -953.4   Mean   :     0.0  
 3rd Qu.: 94194   3rd Qu.: 89084   3rd Qu.:  2952.5   3rd Qu.:  6619.6  
 Max.   :139838   Max.   :101417   Max.   : 33442.4   Max.   :  9648.0

# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que los valores extremos tienen gran influencia y se pronosticará un aumento del consumo.

Tercera iteracción: Si añado \(\gamma\) podré ver la estacionalidad

tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                                alpha=TRUE,
                                beta=TRUE,
                                gamma=TRUE
                                )
summary(tsSM1_HW_Semanal$fitted)

      xhat            level            trend              season        
 Min.   : 36880   Min.   : 67975   Min.   :-32119.9   Min.   :-19389.3  
 1st Qu.: 73545   1st Qu.: 80568   1st Qu.: -6256.7   1st Qu.: -1906.8  
 Median : 86848   Median : 87483   Median :  -534.6   Median :   336.8  
 Mean   : 83905   Mean   : 84859   Mean   :  -953.4   Mean   :     0.0  
 3rd Qu.: 94194   3rd Qu.: 89084   3rd Qu.:  2952.5   3rd Qu.:  6619.6  
 Max.   :139838   Max.   :101417   Max.   : 33442.4   Max.   :  9648.0

# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 7 observaciones, es decir, la predicción del consumo energético semanal total para el año 2010.

## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=7)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha ",
     main="Predicción del consumo energético semanal de la cocina \n Año 2010" )

Vemos como hay un pronóstico de consumo energético mucho mayor que para los años anteriores.

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=7, level = c(95))
## Plot only the forecasted area

plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético semanal para el año 2010", start(2010))

Lavadero

tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,1]
plot.ts(tsSM2_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Evolución del consumo energético según día de la semana \n del lavadero',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.71849, p-value = 4.722e-05

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE,t.window = 7))

Tendencia decreciente del consumo.

Suavizado exponencial de HoltWinters

Primera interacción: rango de predicción

tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                               alpha = TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)

      xhat            level            trend           season      
 Min.   :  5231   Min.   : 17137   Min.   :-3959   Min.   :-20239  
 1st Qu.: 65759   1st Qu.: 80171   1st Qu.:-3959   1st Qu.:-13862  
 Median : 95066   Median : 98804   Median :-3959   Median : -7947  
 Mean   : 89219   Mean   : 93179   Mean   :-3959   Mean   :     0  
 3rd Qu.:109099   3rd Qu.:108252   3rd Qu.:-3959   3rd Qu.:  5329  
 Max.   :161540   Max.   :146103   Max.   :-3959   Max.   : 51483

plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronóstico de consumo energético para el lavadero es de 5231 a 161540 vatios hora.

Segunda interacción: tendencia

tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                               alpha = TRUE, beta=TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)

      xhat            level            trend            season      
 Min.   :-67191   Min.   : 17137   Min.   :-76382   Min.   :-20239  
 1st Qu.: 46662   1st Qu.: 80171   1st Qu.:-20120   1st Qu.:-13862  
 Median : 89948   Median : 98804   Median :-10241   Median : -7947  
 Mean   : 91445   Mean   : 93179   Mean   : -1734   Mean   :     0  
 3rd Qu.:142413   3rd Qu.:108252   3rd Qu.: 11518   3rd Qu.:  5329  
 Max.   :262749   Max.   :146103   Max.   :128966   Max.   : 51483

plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Serie muy afectada por valores extremos, pero se aprecia una tendencia decreciente del consumo energético.

Tercera interacción

tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                               alpha = TRUE,
                               beta=TRUE , gamma=TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)

      xhat            level            trend            season      
 Min.   :-67191   Min.   : 17137   Min.   :-76382   Min.   :-20239  
 1st Qu.: 46662   1st Qu.: 80171   1st Qu.:-20120   1st Qu.:-13862  
 Median : 89948   Median : 98804   Median :-10241   Median : -7947  
 Mean   : 91445   Mean   : 93179   Mean   : -1734   Mean   :     0  
 3rd Qu.:142413   3rd Qu.:108252   3rd Qu.: 11518   3rd Qu.:  5329  
 Max.   :262749   Max.   :146103   Max.   :128966   Max.   : 51483

plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                             #  alpha = TRUE,
                               beta=TRUE , gamma=TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)

      xhat            level            trend            season        
 Min.   : 39417   Min.   : 60834   Min.   :-11558   Min.   :-23890.7  
 1st Qu.: 58672   1st Qu.: 72658   1st Qu.: -7955   1st Qu.:-14149.9  
 Median : 88846   Median : 85738   Median : -6115   Median : -7378.4  
 Mean   : 85244   Mean   : 90171   Mean   : -3035   Mean   : -1891.6  
 3rd Qu.:101345   3rd Qu.:106734   3rd Qu.:  3060   3rd Qu.:  -584.4  
 Max.   :169618   Max.   :128927   Max.   : 10093   Max.   : 51482.6

plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que esta configuración se ajuste mejor.

Pronóstico de HoltWinters para el año 2010

Vamos a predecir el consumo energético de los próximos 7 días, es decir, el consumo energético de cada día de la semana del año 2010

## HoltWinters forecast & plot
     tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=7)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha",
     main="Predicción del consumo energético del lavadero \ según día de la semana del año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=7, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha ", main="Predicción semanal para el año 2010", start(2010))

AC y termo eléctrico

tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,1]
plot.ts(tsSM3_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Evolución según día de la semana del consumo energético del AC y termo \n Serie sin componente estacional',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.5979, p-value = 1.871e-06

Los resíduos no siguen una distribución normal. La descomposición no es buena. Vamos a volver a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE,t.window = 7))

Suavizado exponencial de HoltWinters

Primera interacción: rango de pronóstico

tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                             alpha = TRUE)
summary(tsSM3_HW_Semanal$fitted)

      xhat            level            trend          season        
 Min.   :430345   Min.   :427107   Min.   :3805   Min.   :-14387.8  
 1st Qu.:453104   1st Qu.:446399   1st Qu.:3805   1st Qu.:-10910.1  
 Median :467179   Median :468759   Median :3805   Median :  -566.8  
 Mean   :480676   Mean   :476871   Mean   :3805   Mean   :     0.0  
 3rd Qu.:494414   3rd Qu.:500565   3rd Qu.:3805   3rd Qu.:  8632.7  
 Max.   :566292   Max.   :545292   Max.   :3805   Max.   : 17195.5

plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Predicción del consumo para el termo eléctrico y aire acondicionado entre 430345 y 566292 vatios hora.

Segunda interacción: tendencia de consumo

tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                             alpha = TRUE,beta=TRUE)
summary(tsSM3_HW_Semanal$fitted)

      xhat            level            trend            season        
 Min.   :419999   Min.   :427107   Min.   :-34437   Min.   :-14387.8  
 1st Qu.:447423   1st Qu.:446399   1st Qu.: -7413   1st Qu.:-10910.1  
 Median :462529   Median :468759   Median :  5721   Median :  -566.8  
 Mean   :485585   Mean   :476871   Mean   :  8714   Mean   :     0.0  
 3rd Qu.:507335   3rd Qu.:500565   3rd Qu.: 27545   3rd Qu.:  8632.7  
 Max.   :600575   Max.   :545292   Max.   : 51106   Max.   : 17195.5

plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Se observa una tendencia creciente.

Tercera iteracción

tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                             alpha = TRUE,
                             beta=TRUE, gamma = TRUE)
summary(tsSM3_HW_Semanal$fitted)

      xhat            level            trend            season        
 Min.   :419999   Min.   :427107   Min.   :-34437   Min.   :-14387.8  
 1st Qu.:447423   1st Qu.:446399   1st Qu.: -7413   1st Qu.:-10910.1  
 Median :462529   Median :468759   Median :  5721   Median :  -566.8  
 Mean   :485585   Mean   :476871   Mean   :  8714   Mean   :     0.0  
 3rd Qu.:507335   3rd Qu.:500565   3rd Qu.: 27545   3rd Qu.:  8632.7  
 Max.   :600575   Max.   :545292   Max.   : 51106   Max.   : 17195.5

plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 7 observaciones, es decir, el consumo energético de cada día de la semana para el año 2010.

## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=7)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha",
     main="Predicción del consumo energético semanal \n del AC y termo eléctrico. Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=7, level=c(95))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal del consumo energético para el año 2010", start(2010), ylim=c(0,450000))

Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad

En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.

DatosRealesSeriesSemanales<-Granularidad_semanas  %>% filter(year==2010)

RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_3)[4]


ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranSemanalFor,
               RMSE_SM2_GranSemanalFor,
               RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE
Cocina	28600.80
Lavadero	51998.21
AC y TermoE	218301.40

Los errores son demasiado altos, las predicciones no serán del todo fiables y por tanto no podemos ofrecerlas al cliente.

Análisis de series temporales con metodología Box-Jenkis

Paso 1. Lectura y representación gráfica de los datos

# tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,         frequency=7, start=c(2007,1))
autoplot(tsSM3_Semanal, ts.colour = 'blue', 
         xlab = "Semana del año",
         ylab = "Energía (Varios-hora)", 
         main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

Paso 2. Transformaciones para que la varianza sea estable en el tiempo

library(TSA)
bc=BoxCox.ar(y=tsSM3_Semanal);

bc$mle

[1] 2

bc$ci # El 1 no está en el intervalo, hay que transformar los datos

[1] -1.1  2.0

# El 0 no está en el intervalo, no es transformación logarítmica

Usando la familia de transformaciones BoxCox, nos sugiere un valor de lambda = 2. El valor 1 (no transformar) está en el intervalo, por lo que no hay que hacer una transformación de los datos, podemos asumir que la varianza es constante en el tiempo

# TransftsSM3_Semanal<- tsSM3_Semanal^2

Paso 3. Transformaciones para que la media sea estable en el tiempo

acf(tsSM3_Semanal, main="FAS del consumo de Energía del AC y Termo", lag.max=50)

ndiffs( tsSM3_Semanal) # Necesito una diferencia regular

[1] 1

nsdiffs(tsSM3_Semanal)

[1] 0

Diferencia estacional:

tsSM3_SemanalDif=diff(tsSM3_Semanal, lag=1, diff=1)
acf(tsSM3_SemanalDif, main="FAS de primera diferencia regular del consumo energético", lag=50)

ndiffs(tsSM3_SemanalDif)

[1] 0

Parece la FAS de un Ruido Blanco.

Paso 4. Contrastar la estacionariedad

Test de raiz unitaria.

library(tseries)
adf.test(tsSM3_SemanalDif)


    Augmented Dickey-Fuller Test

data:  tsSM3_SemanalDif
Dickey-Fuller = -6.2783, Lag order = 2, p-value = 0.01
alternative hypothesis: stationary

No existen evidencias significativas para afirmar que los datos no sean estacionarios. La media es estable en el tiempo.

autoplot(tsSM3_SemanalDif, ts.colour = 'blue', 
         xlab = "Semana del año",
         ylab = "Energía (Varios-hora)", 
         main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

Paso 5. Identificar la estructura ARIMA

acf(tsSM3_SemanalDif, main="FAS tras una diferencia regular", lag=50)

Ruido blanco

pacf(tsSM3_SemanalDif, main="FAP tras una diferencia regular", lag=50)

AR(1).

Pasos 6 y 7. Estimación de parámetros y diagnóstico

modelo ARIMA(1,1,0)xARIMA(0,0,0)

ajuste1=arima(tsSM3_Semanal, order=c(1,1,0), seasonal=list(order=c(0,0,0), period=7)); ajuste1


Call:
arima(x = tsSM3_Semanal, order = c(1, 1, 0), seasonal = list(order = c(0, 0, 
    0), period = 7))

Coefficients:
          ar1
      -0.4254
s.e.   0.2312

sigma^2 estimated as 2.967e+09:  log likelihood = -246.59,  aic = 495.17

checkresiduals(ajuste1)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)
Q* = 14.368, df = 3, p-value = 0.002445

Model df: 1.   Total lags used: 4

No pasa la diagnosis: p-value = 0.002445. El modelo no es válido.

ajuste2=auto.arima(tsSM3_Semanal, d = 1,D = 0); ajuste2

Series: tsSM3_Semanal 
ARIMA(3,1,0) with drift 

Coefficients:
          ar1      ar2      ar3     drift
      -1.2341  -1.1548  -0.6579  6092.903
s.e.   0.2171   0.1892   0.2114  1912.826

sigma^2 = 1.26e+09:  log likelihood = -237.16
AIC=484.31   AICc=488.6   BIC=489.29

Modelo propuesto: ARIMA(3,1,0)

aicAutoArima = 496.21 AICModelo1=484.31

checkresiduals(ajuste2)


    Ljung-Box test

data:  Residuals from ARIMA(3,1,0) with drift
Q* = 2.8112, df = 3, p-value = 0.4217

Model df: 4.   Total lags used: 7

El modelo pasa la diagnosis, p-valor=0.422, por tanto es válido. No existen evidencias significativas en contra de que los resíduos sigan un proceso de ruido blanco.

Nos quedamos con el modelo 2.

Paso 8. Predicción de los resultados

pred=forecast(ajuste2,h = 7)# predicción del año 2010
plot(tsSM3_Semanal, xlim=c(2007, 2011), ylim=c(350000,600000))
lines(pred$mean, col="red")

# plot((forecast(ajuste.bueno,h=7))) #muy similar al anterior

Vamos a calcular el error: 4.9418866^{5}, que es mucho mayor que por el otro método.

Pronóstico de consumo para los próximos 12 meses

Lo haremos con el modelo conseguido primero, no con la metodología de Box-Jenkis, que es peor.

Cocina

# Definición de serie temporal

Granularidad_semanas <- Granularidad_horas %>% group_by(year, weekday) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )
seriesSemanalesCompletas <- Granularidad_semanas

## Creamos un objeto TS
 tsSM1_Semanal<- ts(seriesSemanalesCompletas$Sub_metering_1,frequency=7, start=c(2007,1))
 tsSM2_Semanal<- ts(seriesSemanalesCompletas$Sub_metering_2,frequency=7, start=c(2007,1))
 tsSM3_Semanal<- ts(seriesSemanalesCompletas$Sub_metering_3,frequency=7, start=c(2007,1))

# Descomposición de la serie
Componentes_SM1_SemanalSTL<-  tsSM1_Semanal %>%
  stl( t.window=7,s.window="periodic",  robust=TRUE)
# Eliminación de la estacionalidad
tsSM1_Semanal_detrended <- tsSM1_Semanal- Componentes_SM1_SemanalSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalCocinaSemanal <- forecast(tsSM1_HW_Semanal, h=7)
## Plot only the forecasted area
plot(PrediccionFinalCocinaSemanal, 
     ylab= "Vatios-Hora", xlab="Fecha", 
     main="Predicción el consumo energético de la cocina \n según día de para año 2011", start(2010,12), ylim = c(100000,200000))

Las predicciones son las siguientes:

PrediccionFinalCocinaSemanal[["mean"]]

Time Series:
Start = c(2012, 1) 
End = c(2012, 7) 
Frequency = 7 
[1] 173867.6 190256.8 182292.0 184206.4 199466.0 133717.8 125668.2

Según el pronóstico, el día de la semana que más consumo energético total de energía habrá, será el Viernes. Y por el contrario, el día de la semana en que se registrará el menor consumo energético es el Domingo.

Lavadero

# Descomposición de la serie
Componentes_SM2_SemanalSTL<-  
  tsSM2_Semanal %>%
  stl( t.window=7,s.window="periodic",  robust=TRUE)
# Eliminación de la estacionalidad
tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE
 #  alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalLavaderoSemanal <- forecast(tsSM2_HW_Semanal, h=7)
## Plot only the forecasted area
plot(PrediccionFinalLavaderoSemanal, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético del lavadero \n  Según el día de la semana para los próximos 12 meses", start(2010,12), ylim = c(50000,450000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoSemanal[["mean"]]

Time Series:
Start = c(2012, 1) 
End = c(2012, 7) 
Frequency = 7 
[1]  67555.2 134504.6 295918.7 390770.7 420871.3 451580.6 452544.9

Aire acondicionado y termo eléctrico

# Descomposición de la serie
Componentes_SM3_SemanalSTL<-  tsSM3_Semanal %>%
  stl( t.window=7,s.window="periodic",  robust=TRUE)
# Eliminación de la estacionalidad
tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalAcTSemanal <- forecast(tsSM3_HW_Semanal, h=13)
## Plot only the forecasted area
plot(PrediccionFinalAcTSemanal, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético del aire y termo \n para cada día de la semana", start(2010,12), ylim = c(350000,850000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoSemanal[["mean"]]

Time Series:
Start = c(2012, 1) 
End = c(2012, 7) 
Frequency = 7 
[1]  67555.2 134504.6 295918.7 390770.7 420871.3 451580.6 452544.9

Evolución del consumo energético por semana del año.

Serie temporal y visualización

Granularidad_semanas <- Granularidad_horas %>% group_by(year, week) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )

seriesSemanales <- rbind.data.frame(
   filter(Granularidad_semanas, year== 2007  ),
   filter(Granularidad_semanas, year== 2008  ),
   filter(Granularidad_semanas, year== 2009  )
)

## Creamos un objeto TS
   tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1,         frequency=53, start=c(2007,1))
   tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2,         frequency=53, start=c(2007,1))
   tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,         frequency=53, start=c(2007,1))
tsSM_GAP_Semanal<- ts(seriesSemanales$Global_active_power,    frequency=53, start=c(2007,1))

## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la cocina \n Años 2007-2009")

autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la lavandería \n Años 2007-2009")

autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")

autoplot(tsSM_GAP_Semanal, ts.colour = 'green', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía \n Años 2007-2009")

AC y Termo: podría apreciarse cierta estacionalidad Energía total: muelle

Predicción del año 2010 antes de descomponer la serie

Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.

SM1. Cocina

Modelo

fitSM1_Semanal <- tslm(tsSM1_Semanal ~ trend + season) 
summary(fitSM1_Semanal)


Call:
tslm(formula = tsSM1_Semanal ~ trend + season)

Residuals:
     Min       1Q   Median       3Q      Max 
-10942.7  -1970.9   -182.3   2205.3  12203.0 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10835.442   2693.122   4.023 0.000108 ***
trend          -8.872      8.426  -1.053 0.294745    
season2      6872.206   3753.914   1.831 0.069985 .  
season3      6643.745   3753.942   1.770 0.079662 .  
season4      4273.617   3753.989   1.138 0.257536    
season5      3290.489   3754.056   0.877 0.382751    
season6      4326.029   3754.141   1.152 0.251800    
season7       791.901   3754.245   0.211 0.833347    
season8      -313.227   3754.368  -0.083 0.933669    
season9     -3984.688   3754.509  -1.061 0.290985    
season10     6229.851   3754.670   1.659 0.100054    
season11     2045.390   3754.850   0.545 0.587092    
season12     5559.929   3755.048   1.481 0.141693    
season13     3915.802   3755.266   1.043 0.299459    
season14     1586.341   3755.502   0.422 0.673594    
season15     3100.547   3755.757   0.826 0.410934    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4598 on 105 degrees of freedom
Multiple R-squared:  0.4836,    Adjusted R-squared:  0.223 
F-statistic: 1.855 on 53 and 105 DF,  p-value: 0.003663

Resultados:

Valor del coeficiente de determinación: 0.4836. Es un valor bajo, el modelo únicamente explica el 48% de la variabilidad total.
Hay un coeficiente significativamente no nulo.

Predicción IC

Predicción para los próximos 20 días

forecastfitSM1_Semanal <- forecast(fitSM1_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético semanal total del año 2010 \n  a través de un modelo de regresión")

SM2. Lavadero

Modelo

fitSM2_Semanal <- tslm(tsSM2_Semanal ~ trend + season) 
summary(fitSM2_Semanal)


Call:
tslm(formula = tsSM2_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-12710  -2492    262   2070   9448 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  19974.287   2851.564   7.005 2.45e-10 ***
trend          -46.079      8.921  -5.165 1.15e-06 ***
season2      -4864.254   3974.765  -1.224 0.223773    
season3        714.159   3974.795   0.180 0.857756    
season4       -986.429   3974.845  -0.248 0.804489    
season5      -1733.682   3974.915  -0.436 0.663619    
season6       -145.270   3975.005  -0.037 0.970917    
season7      -4670.524   3975.115  -1.175 0.242676    
season8      -1727.111   3975.245  -0.434 0.664841    
season9      -5223.698   3975.395  -1.314 0.191707    
season10      3920.048   3975.565   0.986 0.326383    
season11      1273.127   3975.756   0.320 0.749436    
season12     -1349.793   3975.966  -0.339 0.734920    
season13      -436.381   3976.196  -0.110 0.912818    
season14     -2924.635   3976.446  -0.735 0.463681    
season15     -2526.222   3976.716  -0.635 0.526645    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4868 on 105 degrees of freedom
Multiple R-squared:  0.5326,    Adjusted R-squared:  0.2966 
F-statistic: 2.257 on 53 and 105 DF,  p-value: 0.0001984

Resultados:

Valor del coeficiente de determinación: 0.5326. Es un valor aceptable, el modelo únicamente explica el 53% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM2_Semanal <- forecast(fitSM2_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

SM3. AC y termo

Modelo

fitSM3_Semanal <- tslm(tsSM3_Semanal ~ trend + season) 
summary(fitSM3_Semanal)


Call:
tslm(formula = tsSM3_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-35933  -5644    -17   6305  26151 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  56824.99    7005.03   8.112 9.94e-13 ***
trend           95.69      21.92   4.366 2.97e-05 ***
season2      13899.65    9764.23   1.424  0.15755    
season3      18434.30    9764.31   1.888  0.06180 .  
season4       9370.94    9764.43   0.960  0.33941    
season5      20968.59    9764.60   2.147  0.03406 *  
season6      14556.57    9764.82   1.491  0.13903    
season7       8915.22    9765.09   0.913  0.36335    
season8      -2728.46    9765.41  -0.279  0.78049    
season9     -25352.82    9765.78  -2.596  0.01078 *  
season10      2945.17    9766.20   0.302  0.76358    
season11      8872.81    9766.67   0.908  0.36571    
season12      8371.46    9767.18   0.857  0.39334    
season13     13322.78    9767.75   1.364  0.17550    
season14      9586.09    9768.36   0.981  0.32868    
season15     -1125.59    9769.03  -0.115  0.90849    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11960 on 105 degrees of freedom
Multiple R-squared:  0.7153,    Adjusted R-squared:  0.5715 
F-statistic: 4.977 on 53 and 105 DF,  p-value: 1.341e-12

Resultados:

Valor del coeficiente de determinación: 0.7153. Es un valor que podría ser aceptable, el modelo únicamente explica el 71% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.

Predicción IC

Predicción para los próximos 12 meses (año 2010)

forecastfitSM3_Semanal <- forecast(fitSM3_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="Predicción del consumo energético del año 2010 \n  a través de un modelo de regresión")

Se prevee una bajada del consumo

Comparación de los coeficientes y resultados de cada modelo

###  DatosRealesSeriesSemanales <- Granularidad_semanas %>% filter(year == 2010 ###  )
###  
###  RMSE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
###    R2_SM1_GranSemanal<-0.975
###  MASE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_1)[12]
###  RMSE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
###    R2_SM2_GranSemanal<-0.7717
###  MASE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_2)[12]
###  RMSE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
###    R2_SM3_GranSemanal<-0.7601
###  MASE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_3)[12]
###  
###  ResEvolSemanal<-
###  cbind.data.frame(
###        Submedidor=c("Cocina","Lavadero","AC y TermoE"),
###        RMSE = c(RMSE_SM1_GranSemanal,
###                 RMSE_SM2_GranSemanal,
###                 RMSE_SM3_GranSemanal),
###       # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMen###  sual),
###        R2 = c(R2_SM1_GranSemanal,
###               R2_SM2_GranSemanal,
###               R2_SM3_GranSemanal))
###  ResEvolSemanal %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = ###  "striped")
###  
###  # MASE: Mean absolute error
###  # RMSE: root mean scuare error. Diferencia entre los valores predichos y ###  los observados.

Forecasting descomponiendo la serie (así ok)

Descomponer: quitar tendencia

Descomposición de la serie y visualización

Submetering 1

 tsSM1_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()

 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend             remainder         
 Min.   :-11380.560   Min.   : 9373.389   Min.   :-17306.164  
 1st Qu.: -1525.527   1st Qu.:10411.060   1st Qu.: -1565.107  
 Median :   605.971   Median :11525.838   Median :    50.446  
 Mean   :     0.000   Mean   :11285.767   Mean   :   161.290  
 3rd Qu.:  2468.518   3rd Qu.:12047.600   3rd Qu.:  1583.123  
 Max.   : 12092.727   Max.   :12804.151   Max.   : 18204.879  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     3994         1637      3148          6386 
   %  62.5         25.6      49.3         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8303  0.9452  0.8244  0.9824  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15

#tsSM1_Semanal %>%
#  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
#  autoplot()

Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM1_Semanal, col = 'gray')
lines(Componentes_SM1_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 2

 tsSM2_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()

 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend            remainder         
 Min.   :-11082.099   Min.   :10972.07   Min.   :-17744.999  
 1st Qu.: -2312.887   1st Qu.:11343.85   1st Qu.: -1509.064  
 Median :   306.932   Median :12813.79   Median :   128.053  
 Mean   :     0.000   Mean   :13181.15   Mean   :   101.019  
 3rd Qu.:  2737.120   3rd Qu.:15219.46   3rd Qu.:  1570.218  
 Max.   :  9565.068   Max.   :15803.51   Max.   : 14301.230  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     5050         3876      3079          7916 
   %  63.8         49.0      38.9         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.7600  0.9452  0.7848  0.9927  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15

#tsSM2_Semanal %>%
#  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
#  autoplot()

Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM2_Semanal, col = 'gray')
lines(Componentes_SM2_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Submetering 3

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

tsSM3_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()

 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder        
 Min.   :-47211.03   Min.   :57753.65   Min.   :-56145.14  
 1st Qu.:-10494.71   1st Qu.:58324.34   1st Qu.: -4929.51  
 Median :  4528.45   Median :63141.11   Median :  -321.52  
 Mean   :     0.00   Mean   :62123.83   Mean   :  -707.60  
 3rd Qu.: 11265.72   3rd Qu.:63648.44   3rd Qu.:  3997.01  
 Max.   : 26621.30   Max.   :68148.68   Max.   : 39973.53  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     21760         5324      8927         23820
   %  91.4         22.4      37.5         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8598  0.9452  0.8486  0.9898  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15

Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM3_Semanal, col = 'gray')
lines(Componentes_SM3_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Tendencia creciente del consumo energético

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Sin estacionalidad

Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes

Evolución mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina

#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)

tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,2]
plot.ts(tsSM1_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.83235, p-value = 3.231e-12

Los resíduos no son solo ruido. Aún hay estacionalidad

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                              #  tsSM1_Mensual_Adjusted, 
  beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 53 observaciones, es decir, el consumo energético semanal de cada submedidor para las 53 semanas del año 2010.

## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=53)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético semanal de la cocina \n Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético semanal para el año 2010", start(2010))

Lavadero

## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)

tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,2]
plot.ts(tsSM2_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.88071, p-value = 5.464e-10

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                                # tsSM2_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 7 observaciones.

## HoltWinters forecast & plot
     tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=53)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción semanal del consumo energético del lavadero \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal para el año 2010", start(2010))

AC y termo eléctrico

# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)

tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,2]
plot.ts(tsSM3_Semanal_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.83537, p-value = 4.311e-12

Los resíduos no siguen una distribución normal. La descomposición no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=53)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético semanal del AC y termo eléctrico \n Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal del consumo energético para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad

En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.

RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_3)[4]


ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranSemanalFor,
               RMSE_SM2_GranSemanalFor,
               RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE
Cocina	64689.54
Lavadero	91145.26
AC y TermoE	545792.05

Evolución mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.

Cocina

tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)

Parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="Predicción del consumo energético de la cocina \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Lavadero

## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)

No arece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII))

TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad

Suavizado exponencial de HoltWinters

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=53)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

AC y termo eléctrico

## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)

No parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))

Suavizado exponencial de HoltWinters

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=53)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))

Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad

RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualForII,
               RMSE_SM2_GranMensualForII,
               RMSE_SM3_GranMensualForII)
   
    
)

ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE
Cocina	18462.43
Lavadero	36358.95
AC y TermoE	433409.63

Comparación de resultados. Predicción antes y despues de eliminar estacionalidad

Comparacion<-cbind.data.frame(
  ResEvolMensual2_SinEstacionalidadConTendencia,
  ResEvolMensual2SinTendenciaConEstacionalidad[,2])

colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE Con Estacionalidad	RMSE Sin estacionalidad
Cocina	14278.84	18462.43
Lavadero	44294.27	36358.95
AC y TermoE	325349.90	433409.63

Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.

#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
#             grobs = list(p1,p2),
#           nrow = 2,
#             top="Sin estacionalidad",bottom="Con estacionalidad"
#             )

Evolución diaria del consumo por meses. Mes de Enero

Serie temporal y visualización

En este caso, crearemos una serie para cada mes, ya que la frecuencia de los datos no será la misma, al tener cada mes un número de días diferente.

par(mfrow=c(2,2))
seriesDiariasEnero <- Granularidad_dias %>% filter(year != 2006 & year != 2010 & month == 1)
## Creamos un objeto TS
   tsSM1_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_1, frequency=31, start=c(2007,1))
   
   tsSM2_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_2, frequency=31, start=c(2007,1))
   
   tsSM3_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_3, frequency=31, start=c(2007,1))

## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_DiariaEnero, ts.colour = 'red', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía de la cocina en el mes de Enero \n Años 2007-2009")

autoplot(tsSM2_DiariaEnero, xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía de la lavandería en el mes de Enero \n Años 2007-2009")

autoplot(tsSM3_DiariaEnero, ts.colour = 'blue', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo energético del aire acondicionado y termo en el mes de Enero \n Años 2007-2009")

Parece que no existen muchos patrones que nos podrían ayudar a hacer predicciones, no serían muy fiables. (gráfica muelle)

Forecasting descomponiendo la serie (así ok)

Descomposición de la serie y visualización

Cocina: Submetering 1

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

 tsSM1_DiariaEnero %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()

 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend             remainder        
 Min.   :-1324.873   Min.   : 871.8631   Min.   :-6121.092  
 1st Qu.: -761.466   1st Qu.:1309.3572   1st Qu.: -491.953  
 Median : -158.424   Median :1491.6212   Median :  220.447  
 Mean   :    0.000   Mean   :1522.7948   Mean   :  550.818  
 3rd Qu.:  337.611   3rd Qu.:1775.0069   3rd Qu.:  732.488  
 Max.   : 5397.352   Max.   :2077.1908   Max.   : 9934.040  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     1099.1        465.6    1224.4        1692.0
   %  65.0         27.5      72.4         100.0 

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8124  0.9452  0.7973  0.9880  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 47 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 5 4
 $ inner: int 1
 $ outer: int 15

La componente irregular tiene bastante peso, al igual que la estacional.

Componentes_SM1_DiariaEneroSTL<-
  tsSM1_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_DiariaEneroSTL, ts.colour = 'blue')

La tendencia del consumo energético de la cocina para los meses de Enero es creciente. Se observa una fuerte componenete estacional.

plot.ts( tsSM1_DiariaEnero, col = 'gray',main="Evolución del consumo energético del lavadero \n  de los meses de Enero")
lines(Componentes_SM1_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Lavadero: Submetering 2

 tsSM2_DiariaEnero %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()

 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder        
 Min.   :-1788.921   Min.   :1491.391   Min.   :-5599.893  
 1st Qu.:-1294.941   1st Qu.:1617.887   1st Qu.: -685.829  
 Median : -782.743   Median :1866.286   Median : -122.576  
 Mean   :    0.000   Mean   :1939.768   Mean   :  373.264  
 3rd Qu.:  718.901   3rd Qu.:2339.704   3rd Qu.: 1313.902  
 Max.   : 4203.925   Max.   :2437.168   Max.   : 6567.217  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     2013.8        721.8    1999.7        3494.0
   %  57.6         20.7      57.2         100.0 

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8067  0.9452  0.7789  0.9886  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 47 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 5 4
 $ inner: int 1
 $ outer: int 15

Componentes_SM2_DiariaEneroSTL<- tsSM2_DiariaEnero %>%
  stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_DiariaEneroSTL, ts.colour = 'blue')

No hay una única tendecia, los primeros días del mes de Enero del año 2007 se aprecia un aumento del consumo energético, al ogual que ocurre desde la segunda mitad del mes de enero del año 2008 hasta finales de 2009. Sin embargo, hay un descenso del consumo desde la segunda mitad del año 2007 hasta la misma fecha del año 2008

plot.ts(tsSM2_DiariaEnero, col = 'gray', main="Evolución del consumo energético \n  del lavadero el mes de Enero")
lines(Componentes_SM2_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Aire acondicionado y termo eléctrico: Submetering 3

tsSM3_DiariaEnero %>%
  stl( t.window=31,
      s.window="periodic", 
      robust=TRUE)  %>% summary()

 Call:
 stl(x = ., s.window = "periodic", t.window = 31, robust = TRUE)

 Time.series components:
    seasonal             trend             remainder        
 Min.   :-3370.606   Min.   : 7519.501   Min.   :-9218.843  
 1st Qu.:-1755.521   1st Qu.:10128.318   1st Qu.:-1110.456  
 Median :   -0.644   Median :10857.506   Median : -400.623  
 Mean   :    0.000   Mean   :10612.677   Mean   : -172.452  
 3rd Qu.: 1285.666   3rd Qu.:11340.287   3rd Qu.: 1214.203  
 Max.   : 4774.143   Max.   :11502.819   Max.   : 9935.006  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     3041         1212      2325          3550 
   %  85.7         34.1      65.5         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8385  0.9452  0.8182  0.9845  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 31 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 4 4
 $ inner: int 1
 $ outer: int 15

La componente estacional tiene mucho peso.

Componentes_SM3_DiariaEneroSTL<- tsSM3_DiariaEnero %>%
  stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_DiariaEneroSTL, ts.colour = 'blue')

Se observa un aumento del consumo energético del lavadero a lo largo de los meses de Enero desde los 8000 vatios hasta los 11000.

plot.ts(tsSM3_DiariaEnero, col = 'gray', main="Evolución del consumo energético del AC y termo \n para el mes de Enero")
lines(Componentes_SM3_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Predicción con Holt-Winters (suavizado exponencial) para el año 2010

Evolución mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina

# LE quito la estacionalidad
tsSM1_DiariaEnero_detrended <- tsSM1_DiariaEnero - Componentes_SM1_DiariaEneroSTL$time.series[,1]
plot.ts(tsSM1_DiariaEnero_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Evolución diaria del consumo energético de la cocina \n En el mes de Enero',
        cex.main = 0.85)

Vamos a estudiar los resíduos:

par(mfrow = c(1,2))
  plot(Componentes_SM1_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM1_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM1_DiariaEneroSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM1_DiariaEneroSTL$time.series[, 3]
W = 0.78105, p-value = 1.895e-10

Los resíduos no son solo ruido. Aún hay estacionalidad. Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE, t.window = 31))

Parece que la componente irregular no sigue ningún patrón.

Suavizado exponencial de HoltWinters

Primera interacción

tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
                                   alpha=TRUE)
summary(tsSM1_HW_DiariaEnero$fitted)

      xhat             level           trend           season        
 Min.   :-8059.6   Min.   :-8138   Min.   :1.099   Min.   :-1783.95  
 1st Qu.:  622.2   1st Qu.: 1337   1st Qu.:1.099   1st Qu.:-1110.06  
 Median : 1780.5   Median : 2239   Median :1.099   Median : -715.70  
 Mean   : 2059.4   Mean   : 2058   Mean   :1.099   Mean   :    0.00  
 3rd Qu.: 3900.9   3rd Qu.: 3276   3rd Qu.:1.099   3rd Qu.:  -41.49  
 Max.   :10897.9   Max.   : 7766   Max.   :1.099   Max.   : 9121.65

# optimiza
plot(tsSM1_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronóstico de consumo variará en el rango: 0-10897 vatios.hora.

Segunda interacción:

tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
                                   alpha=TRUE, beta=TRUE)
summary(tsSM1_HW_DiariaEnero$fitted)

      xhat              level           trend               season        
 Min.   :-15283.3   Min.   :-8138   Min.   :-9059.226   Min.   :-1783.95  
 1st Qu.:   247.6   1st Qu.: 1337   1st Qu.:-1401.153   1st Qu.:-1110.06  
 Median :  1412.9   Median : 2239   Median :   -5.467   Median : -715.70  
 Mean   :  2074.9   Mean   : 2058   Mean   :   16.639   Mean   :    0.00  
 3rd Qu.:  4945.5   3rd Qu.: 3276   3rd Qu.: 1526.117   3rd Qu.:  -41.49  
 Max.   : 18063.3   Max.   : 7766   Max.   :13976.323   Max.   : 9121.65

# optimiza
plot(tsSM1_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia del consumo creciente.

Tercera iteracción

tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
                                   alpha=TRUE,
                                   beta=TRUE, gamma=TRUE)
summary(tsSM1_HW_DiariaEnero$fitted)

      xhat              level           trend               season        
 Min.   :-15283.3   Min.   :-8138   Min.   :-9059.226   Min.   :-1783.95  
 1st Qu.:   247.6   1st Qu.: 1337   1st Qu.:-1401.153   1st Qu.:-1110.06  
 Median :  1412.9   Median : 2239   Median :   -5.467   Median : -715.70  
 Mean   :  2074.9   Mean   : 2058   Mean   :   16.639   Mean   :    0.00  
 3rd Qu.:  4945.5   3rd Qu.: 3276   3rd Qu.: 1526.117   3rd Qu.:  -41.49  
 Max.   : 18063.3   Max.   : 7766   Max.   :13976.323   Max.   : 9121.65

# optimiza
plot(tsSM1_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_DiariaEnero_for <- forecast(tsSM1_HW_DiariaEnero, h=31)
plot(tsSM1_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha",
     main="Predicción del consumo energético de la cocina \n Enero Año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_DiariaEnero_forC <- forecast(tsSM1_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción del consumo energético de la cocina \n  Enero del 2010", start(2010))

Lavadero

tsSM2_DiariaEnero_detrended <- tsSM2_DiariaEnero - Componentes_SM2_DiariaEneroSTL$time.series[,1]
plot.ts(tsSM2_DiariaEnero_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Evolución ',
        cex.main = 0.85)

Estudio de los resíduos.

par(mfrow = c(1,2))
  plot(Componentes_SM2_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM2_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM2_DiariaEneroSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM2_DiariaEneroSTL$time.series[, 3]
W = 0.86799, p-value = 1.342e-07

Los resíduos no siguen una distribución normal. La descomposición no es buena, por lo que volvemos a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM2_DiariaEnero_detrended ,s.window = "periodic",
         robust = TRUE,t.window = 31))

Suavizado exponencial de HoltWinters

Primera interacción

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
                              alpha=TRUE)

summary(tsSM2_HW_DiariaEnero$fitted)

      xhat             level             trend            season       
 Min.   :-5033.6   Min.   :-4350.7   Min.   :-19.53   Min.   :-4239.4  
 1st Qu.:  109.6   1st Qu.:  551.7   1st Qu.:-19.53   1st Qu.:-1085.6  
 Median : 1615.7   Median : 2193.8   Median :-19.53   Median : -253.8  
 Mean   : 2167.6   Mean   : 2187.1   Mean   :-19.53   Mean   :    0.0  
 3rd Qu.: 4042.5   3rd Qu.: 2820.3   3rd Qu.:-19.53   3rd Qu.:  903.3  
 Max.   :11795.2   Max.   : 9030.9   Max.   :-19.53   Max.   : 6231.5

plot(tsSM2_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Rango de pronóstico del consumo energético del lavadero: 0-11795 vatios-hora.

Segunda interacción: tendencia

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
                              alpha=TRUE,beta=TRUE)

summary(tsSM2_HW_DiariaEnero$fitted)

      xhat            level             trend                season       
 Min.   :-11600   Min.   :-4350.7   Min.   :-10691.968   Min.   :-4239.4  
 1st Qu.: -2922   1st Qu.:  551.7   1st Qu.: -2362.484   1st Qu.:-1085.6  
 Median :  1808   Median : 2193.8   Median :    -8.537   Median : -253.8  
 Mean   :  2159   Mean   : 2187.1   Mean   :   -27.645   Mean   :    0.0  
 3rd Qu.:  6071   3rd Qu.: 2820.3   3rd Qu.:  2539.702   3rd Qu.:  903.3  
 Max.   : 19223   Max.   : 9030.9   Max.   : 10338.871   Max.   : 6231.5

plot(tsSM2_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Aumento del consumo

Tercera interacción

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
                              alpha=TRUE,
                              beta=TRUE, gamma=TRUE)

summary(tsSM2_HW_DiariaEnero$fitted)

      xhat            level             trend                season       
 Min.   :-11600   Min.   :-4350.7   Min.   :-10691.968   Min.   :-4239.4  
 1st Qu.: -2922   1st Qu.:  551.7   1st Qu.: -2362.484   1st Qu.:-1085.6  
 Median :  1808   Median : 2193.8   Median :    -8.537   Median : -253.8  
 Mean   :  2159   Mean   : 2187.1   Mean   :   -27.645   Mean   :    0.0  
 3rd Qu.:  6071   3rd Qu.: 2820.3   3rd Qu.:  2539.702   3rd Qu.:  903.3  
 Max.   : 19223   Max.   : 9030.9   Max.   : 10338.871   Max.   : 6231.5

plot(tsSM2_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones, es decir, el consumo energético del mes de enero del año 2010.

## HoltWinters forecast & plot
tsSM2_HW_DiariaEnero_for <- forecast(tsSM2_HW_DiariaEnero, h=31)
plot(tsSM2_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del lavadero \n Enero del año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_DiariaEnero_forC <- forecast(tsSM2_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético del lavadero \n Enero del año 2010", start(2010))

AC y termo eléctrico

tsSM3_DiariaEnero_detrended <- tsSM3_DiariaEnero - Componentes_SM3_DiariaEneroSTL$time.series[,1]
plot.ts(tsSM3_DiariaEnero_detrended, 
        ylab = "Energía (Vatios-hora)", 
        main = 'Variación el consumo energético del mes de Enero',
        cex.main = 0.85)

Vamos a estudiar los resíduos:

par(mfrow = c(1,2))
  plot(Componentes_SM3_DiariaEneroSTL$time.series[,3], col = "blue",
       main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM3_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM3_DiariaEneroSTL$time.series[,3])


    Shapiro-Wilk normality test

data:  Componentes_SM3_DiariaEneroSTL$time.series[, 3]
W = 0.90263, p-value = 3.877e-06

Los resíduos no siguen una distribución normal. La descomposición no es buena, por lo que vamos a volver a descomponer la serie:

## Volvemos a descomponer la serie
plot(stl(tsSM3_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE,t.window = 31))

La componente irregular no presenta ningún patrón. Además, observamos una fuerte tendencia creciente del consumo.

Suavizado exponencial de HoltWinters

Primera interacción: rango del pronóstico de consumo.

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_DiariaEnero <- HoltWinters(tsSM3_DiariaEnero_detrended ,
                             alpha=TRUE)
summary(tsSM3_HW_DiariaEnero$fitted)

      xhat           level             trend            season       
 Min.   :-1411   Min.   : -515.4   Min.   :-29.16   Min.   :-7769.7  
 1st Qu.: 8391   1st Qu.: 8869.1   1st Qu.:-29.16   1st Qu.: -971.9  
 Median :10175   Median :10509.4   Median :-29.16   Median : -245.3  
 Mean   :10326   Mean   :10354.9   Mean   :-29.16   Mean   :    0.0  
 3rd Qu.:12873   3rd Qu.:12369.5   3rd Qu.:-29.16   3rd Qu.:  810.6  
 Max.   :22353   Max.   :19582.7   Max.   :-29.16   Max.   : 4510.7

plot(tsSM3_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronóstico de consumo energético del termo eléctrico y el aire acondicionado variará en torno al siguiente rango: 0-22353 vatios-hora

Pronóstico de HoltWinters para el año 2010

Vamos a predecir las próximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_DiariaEnero_for <- forecast(tsSM3_HW_DiariaEnero, h=31)
plot(tsSM3_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="Predicción del consumo energético del AC y termo eléctrico \n Mes de enero del año 2010" )

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_DiariaEnero_forC <- forecast( object = tsSM3_HW_DiariaEnero_for, h=31)
## Plot only the forecasted area
plot(tsSM3_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el mes de Enero del año 2010", start(2010), ylim=c(0,15000))

Comparación de resultados. Predicción despues de eliminar estacionalidad

seriesDiariasEnero <- Granularidad_dias %>% 
  filter( year == 2010 & month == 1)


RMSE_SM1_Enero<-accuracy(tsSM1_HW_DiariaEnero_forC ,seriesDiariasEnero $Sub_metering_1)[4]
RMSE_SM2_Enero<-accuracy(tsSM2_HW_DiariaEnero_forC  ,seriesDiariasEnero $Sub_metering_2)[4]
RMSE_SM3_Enero<-accuracy(tsSM3_HW_DiariaEnero_forC  ,seriesDiariasEnero $Sub_metering_3)[4]



Error<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_Enero,
               RMSE_SM2_Enero,
               RMSE_SM3_Enero)
   
    
)

Error %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")

Submedidor	RMSE
Cocina	150276.386
Lavadero	36593.987
AC y TermoE	7385.421

El error para la cocina es extremadamente alto, sin embargo para el lavadero podría ser aceptable y para el termo eléctrico y aire acondicionado es bastante bajo, el resultado es bueno.

Pronóstico del consumo para los siguientes meses de Enero

Cocina

# Creamos la serie temporal con los datos completos
seriesDiariasEnero <- Granularidad_dias %>% filter( month == 1)
## Creamos un objeto TS
   tsSM1_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_1, frequency=31, start=c(2007,1))
   
   tsSM2_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_2, frequency=31, start=c(2007,1))
   
   tsSM3_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_3, frequency=31, start=c(2007,1))

# Descomponemos la serie
Componentes_SM1_DiariaEneroSTL<-
  tsSM1_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_DiariaEneroSTL, ts.colour = 'blue')

# Eliminación de la estacionalidad
tsSM1_DiariaEnero_detrended <- tsSM1_DiariaEnero - Componentes_SM1_DiariaEneroSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM1_HW_DiariaEnero <- HoltWinters(tsSM1_DiariaEnero_detrended   , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalCocinaEnero <- forecast(tsSM1_HW_DiariaEnero, h=31)
## Plot only the forecasted area
plot(PrediccionFinalCocinaEnero, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético \n  de la cocina para Enero del 2011", start(2011), ylim = c(0,5000))

Las predicciones son las siguientes:

PrediccionFinalCocinaEnero[["mean"]]

Time Series:
Start = c(2011, 1) 
End = c(2011, 31) 
Frequency = 31 
 [1]   3817.7251   1703.2496   1798.5550   1789.3679   1740.9831   1714.8134
 [7]    191.6482   1149.2603   -742.6860  -1190.2384  -1802.0016   4258.6795
[13]  -2247.3347    793.1571  -3181.4382  -3742.3316  -2512.3870  -3889.7260
[19]  -4244.3764   -788.0276   4649.3493  -4534.0677  -7381.3916  -5918.3862
[25]  -6165.6491  -7019.4282  -6945.7927  -6447.0977  -8353.3125  -9155.1739
[31] -10117.8606

Vimos que estas predicciones tenían un error muy grande por lo que no son muy fiables.

Lavadero

# Descomponemos la serie
Componentes_SM2_DiariaEneroSTL<-
  tsSM2_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_DiariaEneroSTL, ts.colour = 'blue')

# Eliminación de la estacionalidad
tsSM2_DiariaEnero_detrended <- tsSM2_DiariaEnero - Componentes_SM2_DiariaEneroSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended   , 
                                seasonal = "additive",
   beta=TRUE,
 gamma=TRUE
#alpha=TRUE
)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalLavaderoEnero <- forecast(tsSM2_HW_DiariaEnero , h=31)
## Plot only the forecasted area
plot(PrediccionFinalLavaderoEnero, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético \n  del lavadero para Enero del 2011", start(2011), ylim=c(0,10000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoEnero[["mean"]]

Time Series:
Start = c(2011, 1) 
End = c(2011, 31) 
Frequency = 31 
 [1] -3208.6987 -2451.0215   337.0196 -8998.2347 -3419.9241 -2154.1838
 [7] -4699.1157 -4641.6730 -2098.9702 -5915.8668 -5908.5963 -3312.5989
[13] -4924.5015 -5895.6914 -5117.6370 -7699.2853 -9518.6745 -7476.8637
[19] -6785.2376 -6129.8531 -8420.1070 -7299.7906 -3701.8091 -6040.7288
[25] -7146.2346 -5477.4638 -2360.3339 -6951.1231 -4149.3855  2215.0182
[31] -7088.5707

Predicciones negativas.

Aire acondicionado y termo eléctrico

# Descomponemos la serie
Componentes_SM3_DiariaEneroSTL<-
  tsSM3_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_DiariaEneroSTL, ts.colour = 'blue')

# Eliminación de la estacionalidad
tsSM3_DiariaEnero_detrended <- tsSM3_DiariaEnero - Componentes_SM3_DiariaEneroSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM3_HW_DiariaEnero <- HoltWinters(tsSM3_DiariaEnero_detrended   , 
        seasonal = "additive",
 #  beta=TRUE
  #gamma=TRUE
  alpha=TRUE
)
plot(tsSM3_HW_DiariaEnero)

# Pronóstico para los siguientes 13 meses 
PrediccionFinalACTermoEnero <- forecast(tsSM3_HW_DiariaEnero , h=31)
## Plot only the forecasted area
plot(PrediccionFinalACTermoEnero, ylab= "Vatios-Hora", xlab="Fecha", main="Predicción el consumo energético \n  del AC y Termo para Enero del 2011", start(2011), ylim = c(0,20000))

Las predicciones son las siguientes:

PrediccionFinalACTermoEnero[["mean"]]

Time Series:
Start = c(2011, 1) 
End = c(2011, 31) 
Frequency = 31 
 [1] 17286.890  9397.946  7612.869 10384.247 10971.474  8229.004  7210.607
 [8]  8490.469  7574.731  6991.382  8002.352  7186.606  6843.313  6456.335
[15]  9699.005  6319.642 10407.530  9105.477  4273.980  9062.211  7955.532
[22]  9057.087  1206.155 10639.813  7764.096  5200.978 10242.690  7921.854
[29]  9748.997  6640.078 11432.974

Predicciones negativas.